This is a jupyter notebook. The "real" notebook can be found here.
Some imports for plotting
import plotly.offline as plotly
import plotly.graph_objs as go
import plotly.figure_factory as ff
#Necessary for the notebook
plotly.init_notebook_mode(connected=True)
algebra
import numpy as np
stuff for the animation
import ipywidgets as widgets
from ipywidgets import interact
and finallypolyfempy
import polyfempy as pf
This code is not particularly interesting.
It converts a tet-mesh (4 indices) into a face-based tri mesh and uses plotly Mesh3d to plot it.
def create_plot_mesh_and_layout(vertices, connectivity, function):
x = vertices[:,0]
y = vertices[:,1]
if vertices.shape[1] == 3:
z = vertices[:,2]
else:
z = np.zeros(x.shape, dtype=x.dtype)
if connectivity.shape[1] == 3:
f = connectivity
else:
#Convert a tet-mesh into face based triangles.
f = np.ndarray([len(connectivity)*4, 3], dtype=np.int64)
for i in range(len(connectivity)):
f[i*4+0] = np.array([connectivity[i][1], connectivity[i][0], connectivity[i][2]])
f[i*4+1] = np.array([connectivity[i][0], connectivity[i][1], connectivity[i][3]])
f[i*4+2] = np.array([connectivity[i][1], connectivity[i][2], connectivity[i][3]])
f[i*4+3] = np.array([connectivity[i][2], connectivity[i][0], connectivity[i][3]])
mesh = go.Mesh3d(x=x, y=y, z=z,
i=f[:,0], j=f[:,1], k=f[:,2],
intensity=function, flatshading=function is not None,
colorscale='Viridis',
contour=dict(show=True, color='#fff'),
hoverinfo="all")
layout = go.Layout(scene=go.layout.Scene(
aspectmode='data',
xaxis=dict(
autorange=True,
showgrid=False,
zeroline=False,
showline=False,
ticks='',
showticklabels=False
),
yaxis=dict(
autorange=True,
showgrid=False,
zeroline=False,
showline=False,
ticks='',
showticklabels=False
),
zaxis=dict(
autorange=True,
showgrid=False,
zeroline=False,
showline=False,
ticks='',
showticklabels=False
)
))
return mesh, layout
def plot(vertices, connectivity, function, camera=None):
mesh, layout = create_plot_mesh_and_layout(vertices, connectivity, function)
if camera is not None:
layout.scene.camera = camera
fig = go.Figure(data=[mesh], layout=layout)
plotly.iplot(fig)
Creates a quad mesh of n_pts x n_pts in the form of a regular grid
def create_quad_mesh(n_pts):
extend = np.linspace(0,1,n_pts)
x, y = np.meshgrid(extend, extend, sparse=False, indexing='xy')
pts = np.column_stack((x.ravel(), y.ravel()))
faces = np.ndarray([(n_pts-1)**2, 4],dtype=np.int32)
index = 0
for i in range(n_pts-1):
for j in range(n_pts-1):
faces[index, :] = np.array([
j + i * n_pts,
j+1 + i * n_pts,
j+1 + (i+1) * n_pts,
j + (i+1) * n_pts
])
index = index + 1
return pts, faces
Set the mesh path
mesh_path = "plane_hole.obj"
create settings
settings = pf.Settings()
pick linear $P_1$ elements. If the mesh would be a quad it would be $Q_1$
settings.discr_order = 1
normalize the mesh to be in $[0,1]^2$
settings.normalize_mesh = True
and choose Young's modulus and poisson ratio
settings.set_material_params("E", 210000)
settings.set_material_params("nu", 0.3)
We are use a linear material model
settings.tensor_formulation = pf.TensorFormulations.LinearElasticity
Next we setup the problem
problem = pf.GenericTensor()
sideset 1 has zero displacement in $x$
problem.add_dirichlet_value(1, [0, 0], [True, False])
sideset 4 has zero displacement in $y$
problem.add_dirichlet_value(4, [0, 0], [False, True])
sideset 3 has a force (Neumann) of [100, 0] applied
problem.add_neumann_value(3, [100, 0])
fianally set the problem
settings.set_problem(problem)
Solve!
solver = pf.Solver()
solver.settings(str(settings))
solver.load_mesh(mesh_path)
solver.solve()
Get the solution
[pts, tets, disp] = solver.get_sampled_solution()
diplace the mesh
vertices = pts + disp
and get the stresses
mises, _ = solver.get_sampled_mises_avg()
finally plot with the above code
top_camera = dict(
up=dict(x=0, y=1, z=0),
center=dict(x=0, y=0, z=0),
eye=dict(x=0, y=0, z=0.8)
)
plot(vertices, tets, mises, top_camera)
Note that we used get_sampled_mises_avg to get the Von Mises stresses because they depend on the Jacobian of the displacement which, becase we use $P_1$ elements, is piece-wise constant. To avoid that effect in get_sampled_mises_avg the mises are averaged per vertex weighted by the area of the triangles. If you want the "real" mises just call
mises = solver.get_sampled_mises()
plot(vertices, tets, mises, top_camera)
Non-linear example. We want to torque a 3D bar around the $z$ direction.
The example is really similar as the one just above.
Sets the mesh and create the settings. As before
mesh_path = "square_beam_h.HYBRID"
settings = pf.Settings()
It is an hex-mesh so we are using $Q_1$
settings.discr_order = 1
In this case the mesh has already the correct size.
settings.normalize_mesh = False
Choose Young's modulus and poisson ratio, as before
settings.set_material_params("E", 200)
settings.set_material_params("nu", 0.35)
Differently from before we want a non-linear material model: NeoHookean
settings.tensor_formulation = pf.TensorFormulations.NeoHookean
and we want to do 5 steps of incremental loading to avoid ambiguities in the rotation
settings.nl_solver_rhs_steps = 5
Now we setup problem with fixed sideset, rotating an number of tours
problem = pf.Torsion()
sideset 5 is fixed
problem.fixed_boundary = 5
sideset 6 rotates
problem.turning_boundary = 6
around the $z$-axis (2)
problem.axis_coordiante = 2
by half a tour
problem.n_turns = 0.5
Now we choose a coarse visualization mesh
settings.vismesh_rel_area = 0.001
and set the problem and solve
settings.set_problem(problem)
#solving!
solver = pf.Solver()
solver.settings(str(settings))
solver.load_mesh(mesh_path)
solver.solve()
takes approx 1 min because it is a complicated non-linear problem!
Get solution and stesses as before
Since we want to show only the surface there is no need of getting the whole volume, so we set boundary_only to True
[pts, tets, disp] = solver.get_sampled_solution(boundary_only=True)
vertices = pts + disp
mises, _ = solver.get_sampled_mises_avg(boundary_only=True)
and plot the 3d result!
plot(vertices, tets, mises)
Create the mesh using the utility function
pts, faces = create_quad_mesh(50)
create settings
settings = pf.Settings()
settings.vismesh_rel_area = 0.001
pick linear $Q_2$ elements for velocity and $Q_1$ for pressure
settings.discr_order = 2
settings.pressure_discr_order = 1
Set the viscosity of the fluid
settings.set_material_params("viscosity", 1)
We select stokes as material model
settings.tensor_formulation = pf.TensorFormulations.Stokes
The default solver do not support mixed formulation, we need to choose UmfPackLU
settings.set_advanced_option("solver_type", "Eigen::UmfPackLU")
We use the standard Driven Cavity problem
problem = pf.DrivenCavity()
we set the problem
settings.set_problem(problem)
We create the solver and set the settings
solver = pf.Solver()
solver.settings(str(settings))
This time we are using pts and faces instead of loading from the disk
solver.set_mesh(pts, faces)
Solve!
solver.solve()
We now get the solution and the pressure
[pts, tris, velocity] = solver.get_sampled_solution()
and plot it!
top_camera = dict(
up=dict(x=0, y=1, z=0),
center=dict(x=0, y=0, z=0),
eye=dict(x=0, y=0, z=1.2)
)
plot(pts, tris, np.linalg.norm(velocity, axis=1), top_camera)
n_pts = len(pts)
scaling = 3
quiver = ff.create_quiver(
pts[0:n_pts:10,0], pts[0:n_pts:10,1],
scaling*velocity[0:n_pts:10,0], scaling*velocity[0:n_pts:10,1])
quiver.layout.yaxis.scaleanchor="x"
quiver.layout.yaxis.scaleratio=1
plotly.iplot(quiver)
Create the mesh using the utility function
pts, faces = create_quad_mesh(50)
create settings
settings = pf.Settings()
pick linear $Q_1$ elements.
settings.discr_order = 1
and choose Young's modulus and poisson ratio
settings.set_material_params("E", 1)
settings.set_material_params("nu", 0.3)
We are use a linear material model
settings.tensor_formulation = pf.TensorFormulations.LinearElasticity
For efficienty in the browser we select a coarse vis mesh
settings.vismesh_rel_area = 0.001
We simulate from 0 to 10s and 50 steps.
settings.tend = 10
settings.time_steps = 50
Next we setup the problem, this doesnt have any parameters. It will...
problem = pf.Gravity()
we set the problem
settings.set_problem(problem)
We create the solver and set the settings
solver = pf.Solver()
solver.settings(str(settings))
This time we are using pts and faces instead of loading from the disk
solver.set_mesh(pts, faces)
Solve!
solver.solve()
Get the solution and the mises
pts = solver.get_sampled_points_frames()
tris = solver.get_sampled_connectivity_frames()
disp = solver.get_sampled_solution_frames()
mises = solver.get_sampled_mises_avg_frames()
Now we are ready to do the animation
First create a figure widget
top_camera = dict(
up=dict(x=1, y=0, z=0),
center=dict(x=0, y=0, z=0),
eye=dict(x=0, y=0, z=1.2)
)
mesh, layout = create_plot_mesh_and_layout(pts[-1],tris[-1],mises[-1])
mesh.cmin = 0
mesh.cmax = 0.25
layout = go.Layout(scene=go.layout.Scene(
camera=top_camera,
aspectratio = dict( x=1, y=1, z=1 ),
aspectmode = 'manual',
xaxis = dict(range=[-0.1, 1.1]),
yaxis = dict(range=[-0.1, 1])
))
fig = go.FigureWidget(data=[mesh], layout=layout)
then we need the callback
def on_value_change(value):
frame_index = value['new']
fig.data[0].x=pts[frame_index][:,0] + disp[frame_index][:,0]
fig.data[0].y=pts[frame_index][:,1] + disp[frame_index][:,1]
fig.data[0].intensity = mises[frame_index]
finally we create an animation widged
play = widgets.Play(
value=0,
min=0,
max=len(mises)-1,
step=1,
description="Press play",
disabled=False
)
slider = widgets.IntSlider(min=0, max=len(mises)-1)
slider.observe(on_value_change, 'value')
play.observe(play, 'value')
widgets.jslink((play, 'value'), (slider, 'value'))
widgets.VBox((widgets.HBox((play, slider)), fig))